import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
colors = [(0, 0, 0), (1, 0, 0)] # black to red
cm = LinearSegmentedColormap.from_list(
"Custom", colors, N=256)
The complex amplitude of Gaussian beam has the form,
$$ \begin{equation} U(r) = A_{0} \dfrac{W_{0}}{W(z)} exp\left[-\dfrac{\rho^2}{W^2(z)}\right] exp\left[-jkz - jk\dfrac{\rho^2}{2R(z)} + j\zeta(z)\right] \end{equation} $$$ \textbf{Beam width} \; \begin{equation} W(z) = \sqrt{\dfrac{\lambda z_{0}}{\pi}} \sqrt{1 + \left(\dfrac{z}{z_{0}}\right)^2} = W_{0} \sqrt{1 + \left(\dfrac{z}{z_{0}}\right)^2} \\ \end{equation} $ $ \begin{equation} \textbf{Wavefront radius of curvature} \; R(z) = z \left[1 + \left(\dfrac{z_{0}}{z}\right)^2\right] \\ \end{equation} $ $ \begin{equation} \textbf{Gouy phase} \; \zeta(z) = tan^{-1}\left(\dfrac{z}{z_{0}}\right) \\ \end{equation} $ $ \begin{equation} \textbf{Rayleigh range} \; z_{0} = \dfrac{\pi W_{0}^2}{\lambda} \\ \end{equation} $ $ \begin{equation} \textbf{Beam waist} \; W_{0} = \sqrt{\dfrac{\lambda z_{0}}{\pi}} \end{equation} $
wavelength = 632.8e-6 # (mm)
z0 = 1;
W0 = np.sqrt(wavelength * z0 / np.pi)
print(f'Beam waist: {W0:.6f} mm')
print(f'Divergence angle: {2 * np.sqrt(wavelength / (z0 * np.pi)):.6f} rad')
print(f'Depth of focus: {2 * z0:.6f} mm')
x = np.linspace(-3*W0, 3*W0, 1000)
y = np.linspace(-3*W0, 3*W0, 1000)
z = np.linspace(-5*z0, 5*z0, 1000)
xx, yy = np.meshgrid(x, y)
rho = np.sqrt(xx**2 + yy**2)
Beam waist: 0.014192 mm Divergence angle: 0.028385 rad Depth of focus: 2.000000 mm
def W(z, wavelength, z0):
W0 = np.sqrt(wavelength * z0 / np.pi)
return W0 * np.sqrt(1 + (z / z0) ** 2)
_W = W(z, wavelength, z0)
_theta = np.sqrt(wavelength / (z0 * np.pi)) * z
fig = plt.figure(figsize=(8,3))
ax = fig.add_subplot(111)
ax.plot(z, _W, color='blue')
ax.plot(z, -_W, color='blue')
ax.plot(z, _theta, linestyle='--', color='red')
ax.plot(z, -_theta, linestyle='--', color='red')
ax.grid(True)
# ax.set_title(r"Intensity along beam axis ($\rho$=0)", fontsize=20)
ax.set_xlabel("z", fontsize=15)
ax.set_ylabel(r"$W(z)$", fontsize=15)
plt.show()
def R(z, z0):
return z * (1 + (z0 / z) ** 2)
fig = plt.figure(figsize=(8,3))
ax = fig.add_subplot(111)
ax.plot(z, R(z, z0), color='blue')
ax.plot(z, z, linestyle='--', color='red')
ax.grid(True)
# ax.set_title(r"Intensity along beam axis ($\rho$=0)", fontsize=20)
ax.set_xlabel("z", fontsize=15)
ax.set_ylabel(r"$R(z)$", fontsize=15)
ax.set_ylim(-6*z0, 6*z0)
plt.show()
The intensity of the Gaussian beam $I(\rho, z) = \left| U(r) \right|^2$ is,
$$ \begin{align} I(\rho, z) &= \left| A_{0} \right|^2 \left[\dfrac{W_{0}}{W(z)}\right]^2 exp\left[-\dfrac{2\rho^2}{W^2(z)}\right] \\ &= I_{0}\left[\dfrac{W_{0}}{W(z)}\right]^2 exp\left[-\dfrac{2\rho^2}{W^2(z)}\right] \end{align} $$W0 = np.sqrt(wavelength * z0 / np.pi)
W0
0.014192480261642175
def intensity(rho, z, wavelength, z0, I0=1):
W0 = np.sqrt(wavelength * z0 / np.pi)
_W = W(z, wavelength, z0)
return I0 * (W0 / _W) ** 2 * np.exp(-2 * (rho / _W) ** 2)
# return I0 * (W0 / _W) ** 2 * np.exp(-2 * rho ** 2 / _W ** 2)
I1 = intensity(0, z, wavelength, z0, I0=1)
fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(111)
ax.plot(z, I1)
ax.grid(True)
ax.set_title(r"Intensity along beam axis ($\rho$=0)", fontsize=20)
ax.set_xlabel("z", fontsize=15)
ax.set_ylabel(r"$I/I_{0}$", fontsize=15)
plt.show()
Assume z at 0, $z_{0}$, $2z_{0}$
z_list = [0, z0, 2*z0]
fig, ax = plt.subplots(1, 3, figsize=(24,5))
for i, _z in enumerate(z_list):
I = intensity(rho, _z, wavelength, z0, I0=1)
ax[i].imshow(I/I.max(),
cmap=cm,
extent=[x.min(), x.max(), y.min(), y.max()], aspect='equal')
ax[i].set_title(f"z={_z:.3f} mm", fontsize=20)
I = intensity(yy, z, wavelength, z0, I0=1)
fig = plt.figure(figsize=(8,2))
ax = fig.add_subplot(111)
ax.imshow(I, cmap='jet', origin='lower', extent=[z.min(), z.max(), y.min(), y.max()], aspect='auto')
# ax.imshow(I, cmap='jet', origin='lower')
ax.grid(True)
ax.set_title("Intensity at y-z plane", fontsize=20)
ax.set_xlabel("z", fontsize=15)
ax.set_ylabel("y", fontsize=15)
plt.show()
def phase(z, rho, wavelength, z0):
k = 2 * np.pi / wavelength
longi = k * z - np.arctan(z / z0)
trans = k * rho ** 2 / (2 * R(z, z0))
return longi + trans
def amp(z, rho, wavelength, z0, A0=1):
W0 = np.sqrt(wavelength * z0 / np.pi)
_W = W(z, wavelength, z0)
return A0 * (W0 / _W) * np.exp(-(rho / _W) ** 2)
z = np.linspace(-2*z0, 2*z0, 2000)
yy, zz = np.meshgrid(y, z)
U = amp(zz, yy, wavelength, z0, A0=1) * np.exp(-1j * phase(zz, yy, wavelength, z0))
fig = plt.figure(figsize=(12,4))
ax = fig.add_subplot(111)
ax.imshow(np.abs(U.real).T**0.6,
cmap='jet',
origin='lower',
aspect='auto',
extent=[z.min(), z.max(), y.min(), y.max()])
ax.set_xlabel("z", fontsize=15)
ax.set_ylabel(r"$\rho$", fontsize=15)
plt.show()
The complex amplitude of Hermite-Gaussian beam has the form,
$$ \begin{equation} U_{l,m}(x,y,z) = A_{l,m} \left[ \dfrac{W_{0}}{W(z)}\right] \mathbb{G}_{l}\left[\dfrac{\sqrt{2}x}{W(z)}\right] \mathbb{G}_{m}\left[\dfrac{\sqrt{2}y}{W(z)}\right] exp\left[-jkz - jk\dfrac{x^2+y^2}{2R(z)} + j(l+m+1)\zeta(z)\right] \end{equation} $$$ \textbf{Hermite-Gaussian function} \; \begin{equation} \mathbb{G}_{l}(u) = \mathbb{H}_{l}(u) \; exp\left(-\dfrac{u^2}{2}\right), \; l=0,1,2,... \\ \end{equation} $ $ \begin{equation} \textbf{Hermite polynomial} \; \mathbb{H}_{l+1}(u) = 2u\mathbb{H}_{l}(u) - 2l\mathbb{H}_{l-1}(u) ,\:\text{where} \; \mathbb{H}_{0}(u)=1, \; and \; \mathbb{H}_{1}(u)=2u \end{equation} $
def hermite(x, order):
if order == 0:
return np.ones_like(x)
elif order == 1:
return 2 * x
return 2 * x * hermite(x, order-1) - 2 * (order-1) * hermite(x, order-2)
def hermite_gauss(x, order):
h = hermite(x, order)
return h * np.exp(-x ** 2 / 2)
fig, ax = plt.subplots(1, 4, figsize=(32,5))
x = np.linspace(-6, 6, 500)
for i in range(4):
hg = hermite_gauss(x, order=i)
ax[i].plot(x, hg, color='red')
ax[i].fill_between(x, hg, color='#D4D4D4')
ax[i].grid(True)
ax[i].set_title(f'HG with order {i}', fontsize=20)
plt.show()
# for i, _z in enumerate(z_list):
# I = intensity(rho, _z, wavelength, z0, I0=1)
# ax[i].imshow(I/I.max(),
# cmap=cm,
# extent=[x.min(), x.max(), y.min(), y.max()], aspect='equal')
# ax[i].set_title(f"z={_z:.3f} mm", fontsize=20)
wavelength = 632.8e-6 # (mm)
z0 = 1;
W0 = np.sqrt(wavelength * z0 / np.pi)
print(f'Beam waist: {W0:.6f} mm')
print(f'Divergence angle: {2 * np.sqrt(wavelength / (z0 * np.pi)):.6f} rad')
print(f'Depth of focus: {2 * z0:.6f} mm')
x = np.linspace(-3*W0, 3*W0, 1000)
y = np.linspace(-3*W0, 3*W0, 1000)
z = np.linspace(-5*z0, 5*z0, 1000)
xx, yy = np.meshgrid(x, y)
rho = np.sqrt(xx**2 + yy**2)
Beam waist: 0.014192 mm Divergence angle: 0.028385 rad Depth of focus: 2.000000 mm
The intensity of the Hermite-Gaussian beam $I_{l,m} = \left| U_{l,m} \right|^2$ is,
$$ \begin{equation} I_{l,m}(x,y,z) = \left|A_{l,m}\right|^2 \left[ \dfrac{W_{0}}{W(z)}\right]^2 \mathbb{G}_{l}^2\left[\dfrac{\sqrt{2}x}{W(z)}\right] \mathbb{G}_{m}^2\left[\dfrac{\sqrt{2}y}{W(z)}\right] \end{equation} $$def hg_intensity(x, y, z, l, m, wavelength, z0, I0=1):
W0 = np.sqrt(wavelength * z0 / np.pi)
_W = W(z, wavelength, z0)
hg_x = hermite_gauss(np.sqrt(2) * x / _W, l)
hg_y = hermite_gauss(np.sqrt(2) * y / _W, m)
return I0 * (W0 / _W * hg_x * hg_y) ** 2
fig, ax = plt.subplots(4, 4, figsize=(32,32))
for l in range(4):
for m in range(4):
hg = hg_intensity(xx, yy, 0, l, m, wavelength, z0)
ax[m,l].imshow(hg, cmap=cm)
ax[m,l].set_title(f"HG with order {(l,m)}", fontsize=30)
ax[m,l].set_xticks([])
ax[m,l].set_yticks([])
plt.show()
The complex amplitude of Laguerre-Gaussian beam has the form,
$$ \begin{equation} U_{l,m}(\rho,\phi,z) = A_{l,m} \left[ \dfrac{W_{0}}{W(z)}\right] \left( \dfrac{\rho}{W(z)}\right)^{l} \mathbb{L}^{l}_{m}\left(\dfrac{2\rho^{2}}{W^{2}(z)}\right) exp\left(-\dfrac{\rho^{2}}{W^{2}(z)}\right) exp\left[-jkz - jk\dfrac{\rho^{2}}{2R(z)} \mp{jl\phi}+ j(l+2m+1)\zeta(z)\right] \end{equation} $$$ \textbf{Generalized Laguerre polynomial} \; \begin{equation} \mathbb{L}^{l}_{m+1}(u) = \dfrac{(2m+1+l-u)\mathbb{L}^{l}_{m}(u) - (m+l)\mathbb{L}^{l}_{m-1}(u)}{m+1} , \\ \:\text{where} \; \mathbb{L}^{l}_{0}(u) = 1 \; and \; \mathbb{L}^{l}_{1}(u) = 1 + l - u \end{equation} $
def laguerre(x, l, m):
if m == 0:
return np.ones_like(x)
elif m == 1:
return 1 + l - x
else:
return ((2*m-1+l-x)*laguerre(x, l, m-1) - (m-1-l)*laguerre(x, l, m-2)) / m
x = np.linspace(-2, 10, 500)
fig = plt.figure(figsize=(8,6))
ax = plt.subplot(111)
for l in range(4):
for m in range(4):
ax.plot(x, laguerre(x, l, m))
ax.grid(True)
ax.set_ylim(-5, 10)
plt.show()
x = np.linspace(-5, 20, 500)
fig = plt.figure(figsize=(8,6))
ax = plt.subplot(111)
for m in range(6):
ax.plot(x, laguerre(x, 0, m))
ax.grid(True)
ax.set_ylim(-10, 20)
plt.show()
The intensity of the Laguerre-Gaussian beam $I_{l,m} = \left| U_{l,m} \right|^2$ is,
$$ \begin{equation} I_{l,m}(\rho,\phi,z) = \left|A_{l,m}\right|^{2} \left[ \dfrac{W_{0}}{W(z)}\right]^{2} \left( \dfrac{\rho}{W(z)}\right)^{2l} {\mathbb{L}^{l}_{m}}^{2}\left(\dfrac{2\rho^{2}}{W^{2}(z)}\right) exp\left(-\dfrac{2\rho^{2}}{W^{2}(z)}\right) \end{equation} $$def lg_intensity(rho, z, l, m, wavelength, z0, I0=1):
W0 = np.sqrt(wavelength * z0 / np.pi)
_W = W(z, wavelength, z0)
lg = laguerre(2 * (rho / _W) ** 2, l, m)
return I0 * (W0 / _W * lg) ** 2 * np.exp(-2 * (rho / _W) ** 2) * (rho / _W) ** (2*l)
wavelength = 632.8e-6 # (mm)
z0 = 1;
W0 = np.sqrt(wavelength * z0 / np.pi)
print(f'Beam waist: {W0:.6f} mm')
print(f'Divergence angle: {2 * np.sqrt(wavelength / (z0 * np.pi)):.6f} rad')
print(f'Depth of focus: {2 * z0:.6f} mm')
x = np.linspace(-3*W0, 3*W0, 1000)
y = np.linspace(-3*W0, 3*W0, 1000)
z = np.linspace(-5*z0, 5*z0, 1000)
xx, yy = np.meshgrid(x, y)
rho = np.sqrt(xx**2 + yy**2)
Beam waist: 0.014192 mm Divergence angle: 0.028385 rad Depth of focus: 2.000000 mm
fig, ax = plt.subplots(4, 4, figsize=(32,32))
for l in range(4):
for m in range(4):
lg = lg_intensity(rho, 0, l, m, wavelength, z0)
ax[m,l].imshow(lg, cmap='jet')
ax[m,l].set_title(f"LG with order {(l,m)}", fontsize=30)
ax[m,l].set_xticks([])
ax[m,l].set_yticks([])
plt.show()
The complex envelope of Bessel beam has the form,
$$ \begin{equation} A(x,y) = A_{m} J_{m}(k_{T}\rho) e^{jm\phi}, \; \text{where} \; k_{T}^2 + \beta^2 = k^2 \end{equation} $$$ \textbf{Bessel function of the first kind and $m^{th}$ order} \; \begin{equation} J_{m}(x) = \sum_{n=0}^{\infty}{\dfrac{(-1)^{n}}{n! \Gamma(n+m+1)} \left( \dfrac{x}{2}\right)^{2n+m}} \end{equation} $
Thus, the complex amplitude is,
$$ \begin{equation} U(r) = A_{m} J_{m}(k_{T}\rho) e^{jm\phi}, \; \text{where} \; k_{T}^2 + \beta^2 = k^2 \end{equation} $$If $m=0$, the wave has complex amplitude, $ \begin{equation} U(r) = A_{0} J_{0}(k_{T}\rho) e^{-j \beta z} \end{equation} $ ,and the intensity distribution $I(\rho, \phi, z) = \left|A_{0}\right|^{2} J_{0}^{2}(k_{T}\rho)$
import plotly.graph_objects as go
fig= go.Figure(data=go.Isosurface(
x=[0,0,0,0,1,1,1,1],
y=[1,0,1,0,1,0,1,0],
z=[1,1,0,0,1,1,0,0],
value=[1,2,3,4,5,6,7,8],
isomin=2,
isomax=6,
))
fig.show()